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The aim of this paper is two-fold. First, we present a phase space perspective on long range double 
ionization in a one dimensional model of helium. The dynamics is simulated with the periodic von 
Neumann (PvB) method, an exact quantum method based on a lattice of phase space Gaussians. 
Second, we benchmark the method by comparing to the Multiconfiguration Time-dependent Hartree 
method. The PvB approach is found to be faster than the standard MCTDH code for the dynamics 
calculations and to give better accuracy control. 


I. INTRODUCTION 

The possibility of controlling electron dynamics in iso¬ 
lated atoms and molecules has generated growing interest 
in recent years 0,i- Since the timescale for electron dy¬ 
namics is attoseconds to femtoseconds, to achieve control 
one would ideally like to have intense fields (10^"^ — 10^^ 
W/cm^) of several attosecond to subfemtosecond dura¬ 
tion. It is now routinely possible to produce NIR fields of 
10 fs duration with intensity on the order of 10^^ — 10^® 
W/cm^ and XUV pulses of femtosecond duration with 
intensities on the order of 10^^ W/cm^. The desired com¬ 
bination of XUV pulses of sub-fs duration with intensities 
of 10^^ — 10^® W/cm^ is still not readily achievable and 
therefore a slew of recent experiments have employed NIR 
-I- XUV pulses, so that the short and weak XUV pulse is 
boosted by the intensity of the NIR pulse. 

From the point of view of simulation, the situation is 
reversed. The relatively weak XUV pulse, coupling just 
a few angular momentum states, is much easier to sim¬ 
ulate, using e.g. hyperspherical based methods 0,i or 
the R-matrix method [^. In contrast, the intense NIR 
fields can couple scores or hundreds of angular momen¬ 
tum states, making the simulation extremely expensive 
or impossible with hyperspherical based methods unless 
the held intensity in the simulation is reduced well below 
the usual experimental range of values 0| . 

Because of the challenge of hyperspherical based sim¬ 
ulations, several alternative approaches have been ex¬ 
plored. Among them, the Multi-Conhguration Time De¬ 
pendent Hartree method (MCTDH) is especially 

appealing due to its favorable scaling properties. Origi¬ 
nally developed to simulate the vibrational dynamics of 
polyatomic molecules, it has been extended to electronic 
dynamics as well. It was shown that it can describe the 
double ionization in a simple model of helium where each 
electron is described by only one degree of freedom [T7I| . 
More recently, it was extended to describe double ioniza¬ 
tion in a two dimensional and a three dimensional 

model of helium [l^ with NIR excitation and two active 
electrons. In contrast with the hyperspherical approach. 
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MCTDH is able to simulate both strong NIR and XUV 
pulses. 

Several other new methods are currently being de¬ 
veloped to describe multielectron dynamics, e.g. the 
time-dependent generalized-active-^ace configuration- 
interaction (TD-GASCI) method [Q and the B-spline 
algebraic diagrammatic construction A recent re¬ 

view of this field with a focus on MCTDH can be found 



The above-mentioned calculations have given much 
insight into the participation of multiple electrons in 
the ionization process in the presence of strong fields. 
Despite the significant advances in our understanding, 
there are several motivations for new computational ap¬ 
proaches. 1. Certain aspects of the ionization process 
may require a much larger spatial range than the one 
used in the above studies. For example, the study of the 
time delay in [s^ required a grid more than an order of 
magnitude larger than the one used in [T^, [HI, and as 
a result was restricted to the single-active-electron ap¬ 
proximation. Note that the t-SURFF method [s^ can 
describe long range dynamics by construction, but it ne¬ 
glects electron-electron interaction at large distance. The 
phase space approach presented here, though applied to 
two electrons in 1-d, is designed to be more efficient than 
current methods when extremely large grids are required, 
without neglecting any interaction. 2. The classical 
picture underlying high harmonic generation by strong 
NIR pulses is the so-called three-step model consisting of 
strong field tunneling ionization, quasi-free electron prop¬ 
agation and recollision . The quantum analog of this 
three-step model was formulated shortly after the clas¬ 
sical model and captures most of the key features [Hl| . 
Normally, full quantum simulation methods do not ex¬ 
ploit this underlying classical structure. In contrast, the 
phase-space propagation method presented here is a fully 
quantum method that still allows one to see the under¬ 
lying classical structure of the dynamics and therefore to 
understand the predictive limits of the classical model. 3. 
Usually calculations are performed in the coordinate rep¬ 
resentation despite the fact that some of the key experi¬ 
mental observables are associated with the momentum of 
the ionizing electrons. An approach that captures the key 
dynamics simultaneously in coordinate and momentum, 
i.e. in phase space, could have significant advantages in 
the analysis of the correlated wavefunction. The advan- 
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tages of points 2. and 3. could in principle be obtained 
with other phase space representations, e.g. the Wigner 
function, but those representations are generally not as 
compact as the one presented here. 

The present approach, called PvB, is based on a pe¬ 
riodic lattice of phase space Gaussians, called the von 
Neumann basis [l^. The fact that the method employs 
a phase space representation enables one to calculate the 
time evolving wavefunction only in those regions of phase 
space where there is significant wavefunction amplitude. 
A by-product of the method is that it gives intuition into 
the electron dynamics in phase space and the underly¬ 
ing classical correspondence. The PvB method has been 
applied successfully to the ionization process of one di¬ 
mensional hydrogen under a combination of strong NIR 
and XUV pulses [l^. Here we extend the method to 
study the double ionization of two electrons in 1-d in the 
presence of NIR and XUV pulses. In addition to the mul- 
tidimentional formulation, a key feature of the method 
presented here is the algorithm to dynamically adapt the 
phase space during the quantum dynamics. 

The aim of this paper is twofold: first, to present 
the advantages of a phase space point of view for ana¬ 
lyzing double ionization and second, to benchmark the 
PvB method by comparing it with the MCTDH method. 
While the MCTDH algorithm is faster to calculate the 
initial state, for the dynamics calculations presented here, 
the PvB approach is found to be significantly faster than 
the standard Heidelberg MCTDH code, as well as to give 
better accuracy control. However, it is important to men¬ 
tion that this version of MCTDH has not been optimized 
for ionizing systems and that other implementations of 
MCTDH may be significantly more efficient. In the fu¬ 
ture it might be interesting to explore combining the PvB 
approach with MCTDH to obtain additional computa¬ 
tional efficiency. 

The paper is organized as follows. In Section 2 we re¬ 
view the periodic von Neumann (PvB) and the MCTDH 
approaches to simulate quantum dynamics. In Section 
3.1 we present the model system which will be used to 
compare the two approaches. Next we detail the results: 
a comparison of the performance of the two approaches 
for eigen-decomposition (Section 3.2) and ionization dy¬ 
namics, (Section 3.3). Section 3.3 includes a series of 
phase space pictures of double ionization to highlight the 
insight given by this approach. Section 4 concludes with 
a review of the advantages and current limitations of the 
PvB method and an outlook of how it may be further 
developed. 


II. METHODOLOGY 

A. The bi-orthogonal von Neumann Basis 

Before going into details we provide a brief summary 
of the method. The von Neumann (vN) lattice is a lat¬ 
tice of Gaussians in phase space with one Gaussian per 


phase space cell of area h‘^ {h is Planck’s constant and d 
is the dimension) [d^. This basis is known to be com¬ 
plete but not overcomplete on the infinite plane [4l| [? 
]. The appeal of this basis is that the Gaussians can 
be placed only where needed in phase space, and hence 
classical intuition can be used to guide and to interpret 
the quantum calculation. However, in any calculation 
on a truncated phase space the vN lattice is known to 
have severe problems with convergence [H, |4^ . We have 
shown that by modifying the Gaussians gi to be band- 
limited and periodic, one can obtain Fourier accuracy 
[dj . In order to prune basis functions from regions 
of phase space that are not necessary one actually needs 
to use the biorthogonal partners bi instead of the gi as 
the basis functions as discussed below (for a fuller ex¬ 
planation see [alii, [13). We call this “biorthogonal 
exchange”, giving rise to the name of the method ” peri¬ 
odic von Neumann basis with biorthogonal exchange” or 
PvB. An alternative method to converge the truncated 
von Neumann lattice has been developed by Poirier (d^ - 


1. The Underlying Hilbert Space 

We begin by defining the Hilbert space which serves 
as the foundation for all further discussion — a trun¬ 
cated discrete Fourier basis. As we will see, this finite¬ 
dimensional Hilbert space corresponds to a rectangular 
area in phase-space. 

Choosing a finite region of length L in x, we may as¬ 
sume, without loss of generality, cyclic boundary condi¬ 
tions / (x) = f {x + L). This, in turn, implies that the 
functions in this interval are spanned by the orthogonal 
functions of the form exp (2Tri^ri) = exp (ifc„x), Vn £ 
Z, with kn = Next, we limit ourselves to |A:„| < 

ATiim, i.e. n G [-^§- + 1; with K . end¬ 

ing up with a rectangular area of phase-space spanned 
by a discrete number, N = 2-|^, of complex exponential 
functions — a set we shall denote as the spectral basis. 
This Fourier grid of N points spans an area oi Nh in 
phase-space (L x 2P = = Nh). 

Nyquist’s theorem ensures that by sampling the in¬ 
terval at N uniformly-space points, i.e. at resolution 
Sx = ^/n, we can fully reconstruct any function resid¬ 
ing inside the phase-space, i.e. the functions spanned by 
the spectral basis. Given any two of L, K and Sx (or 
N), we may define the set of sampling points, henceforth 
denoted the Fourier grid (FG). 

The set of functions within the Hilbert space that take 
on the value 0 at all grid points except for a single grid 
point where the value is 1, is known as the pseudo-spectral 
basis. These functions can be shown to be orthonormal 
and span the same Hilbert space as the spectral basis. 
For the Fourier grid, the pseudo-spectral functions are 
the N periodic sine functions sinc„ (x) = ; 

centered at the N grid points M- 
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2. PvB - Bi-orthogonal von Neumann Basis 

The Fourier basis described above is capable of giving 
high accuracy but does not generate a sparse represen¬ 
tation of the state. However, given a Fourier basis of 
N functions, corresponding to a rectangular region of 
phase-space, we may span the same space using any set 
of N linearly-independent functions which are themselves 
within that subspace. Specifically, to create this linear 
combination we use a set of phase space Gaussians, cen¬ 
tered on a grid of \/aN x ^^N/a points within the phase 
space spanned by the Fourier grid. See fig. [1] 



extremely delocalized. 

In matrix notation, the biorthogonality takes the form 
= G^H = J, (2) 

i.e. B = Gf\ We may introduce a vector notation 
■05 = = G' 1 ' 0 , where the elements = (fffcl 0) 

are the coefficients of 0 in the B basis. We then find that 
0 = Bi/jb = B (G'^0) = (HG'I') 0, consistent with eq. [2] 


3. The Schrodinger Equation 

Let us derive the form of the Schrodinger equation in 
the bi-orthogonal von Neumann basis. Starting with the 
standard expressions. Hip = i50 and 9*0 = sub¬ 

stituting 0 = BipB and shifting the time-independent B 
to the other side we end up with 

B-^HBPjb = Xe^b . (3) 

Similarly, the time-dependent Schrodinger equation be¬ 
comes 

dtipB = - ^B~^HBipB (4) 

Note the appearance of a similarity transform of the 
Hamiltonian, as opposed to the standard unitary rota¬ 
tion. Also, recall B~^ = G0 


FIG. 1. A 5 X 5 von Neumann grid of Gaussians (real compo¬ 
nent), spanning the same subspace as a 25-point Fourier grid, 
serving as the set of bras for the bi-orthogonal von Neumann 
basis. The figure depicts a mixed representation showing the 
oscillations in the basis functions in a;-space centered at dif¬ 
ferent fcs. 
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f X — X. 

V 2CTa; 




( 1 ) 

We shall term this grid the von Neumann lattice and 
the basis formed by this set of N Gaussians states as the 
von Neumann basis, with each Gaussian loosely speaking 
spanning an area of h. 

Let us denote the set of Gaussian states of the von 
Neumann basis as (/ = {gk }and their matrix repre¬ 
sentation in the Fourier grid as the NxN matrix G, with 
Gj^k = gk (xj), where k is an enumeration of the N Gaus¬ 
sians and Xj are the Fourier Grid points. As the Gaus¬ 
sians are non-orthogonal, we define the bi-orthogonal ba¬ 
sis to G, which we shall term B: B = {bk}^^i where 
{gk\ bj) = Skj. We denote the matrix representation of 
B as the NxN matrix B. Note that while G is com¬ 
posed of well-localized Gaussians, the functions in B are 


4- The Reduced PvB Basis 

As discussed, in the bi-orthogonal von Neumann basis 
the representation of localized wavefunctions is sparse. 
The coefficients that are close to zero can be neglected 
with minimal and well-controlled loss of accuracy and 
therefore this representation can be used to propagate 
dynamics in an efficient manner. Specifically, as |0) = 
J2k i’Bk\bk) with 0B^ = {gk\ 0), we expect the ipB vector 
to have a negligible value for areas in phase-space where 
0 is not present. One may therefore reduce the B ba¬ 
sis to the subset of &-vectors whose coefficients are above 
some wavefunction amplitude cutoff {W). Formally, we 
are projecting the state to a subset of the B basis vectors. 
We shall term this the reduced basis. In matrix represen¬ 
tation, this corresponds to the selection of some columns 
of B, giving a NxR matrix, with R<N, that we denote 
by B. The state is initially represented as a column vec¬ 
tor containing N coefficients. After this reduction the 
state is represented by a column vector of R coefficients 
such that all of them have an absolute value above W. 
This vector is denoted 0^ and termed the reduced state. 
Note that R can be several order of magnitude smaller 
than N. 

As B is no longer square, B~^ is not well defined, and 
it must be replaced in eqs. I3I4I by the pseudo-inverse. 
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defined by B^^ B = I. This gives 

(^B^bY' (^B^HB)i;^ = XE^B ( 5 ) 
dt^B = -\ {b^Y ^B- ( 6 ) 

The truncation of the B basis defines a projection that 
possesses interesting properties detailed in [^. In par¬ 
ticular, it was shown in |l7| that this is the optimal way 
to project a state in a subspace to obtain a reduced state 
as close as possible to the original state. The process pre¬ 
serves the form of the {bk} functions, i.e. bk = bk, but 
replaces the Gaussians {gk} by a new set of deformed 
Gaussians, {gt}, that is orthogonal to the reduced set of 

{h}. 

To solve eq. [S] numerically one can use any standard 
quantum propagation algorithm able to deal with non- 
Hermitian Hamiltonians. For this study we use the Tay¬ 
lor propagation presented in [d^, which we found to be 
the most efficient for our method. 

Note that in this formulation, memory and time re¬ 
quirements scale with the reduced basis size, not the full 
Hilbert space size. 

5. Finding Eigenmodes in PvB 

We now turn to computing the eigenmodes of a given 
potential in the reduced PvB basis. We are faced with 
a seemingly intractable problem. On the one hand, we 
are unable to represent the full Hilbert space in memory 
due to its unmanageable size, and on the other hand, 
we have no foreknowledge of which &-vectors make up 
the appropriate reduced basis. We therefore utilize an 
iterative algorithm, detailed in [2^. The iteration starts 
at the locations most likely to be included in the reduced 
basis — the potential’s local minima — and expands the 
reduced space as needed. We shall refer to regions on the 
von Neumann lattice at which Gaussians are centered as 
cells, allowing us to speak of “neighboring cells in phase- 
space” , “cells at the boundary of the reduced basis set”, 
etc. 

Following is a pseudocode for the iterative algorithm: 

10 Define the initial reduced basis consisting 
of the cells at the position X of the 
potential’s local minima, with p = 0. 

20 Compute eigenmodes for the current reduced 
basis. 

30 If everywhere on the boundary of the 
reduced basis, all eigenmode aunplitudes are 
below the specified accuracy threshold, stop. 

If not, continue. 

40 Remove all phase-space cells where the 
cunplitude is below the wavefunction accuracy 


cutoff. 

50 Expand the reduced basis to all neighboring 
cells (i.e. all cells at or below some 
distcince r from the current reduced basis) 

60 Compute the additional elements of the now 
expanded, reduced Hcuniltonian. 

70 Go to 20 


A more detailed description of this, and the following 
dynamics algorithm, is given in [^ . 

6. Dynamics in the Reduced PvB Basis 

The main difference between the dynamics algorithm 
used in this work and the one used in [l^ , beyond adap¬ 
tation to multidimensional systems, is a new dynamic 
method for choosing the reduced basis set. In [1^, the 
total time was divided in 8 segments. For each of these 
segments some concatenation of rectangles of phase space 
was chosen as the reduced basis set based on classical tra¬ 
jectories. Moreover, the dynamics was computed with a 
fixed time step. Here, the methodology is refined based 
on the following insight: given that the time evolution of 
the wavefunction is continuous in the phase-space, we 
tailor the reduced basis as the wavepacket evolves in 
time. We do this by monitoring the wavefunction am¬ 
plitude at the boundary of the reduced space: If it rises 
above the specified accuracy threshold, for example at 
the “bow” of a travelling wave-packet, the reduced basis 
is expanded in the region of the boundary. Gonversely, 
as the amplitude falls below the threshold at the wave- 
packet’s “stern”, vectors are removed from the reduced 
basis. Finally, we dynamically adapt the time step to 
the speed of the wavepacket at the boundary. If the am¬ 
plitudes at the boundaries grow faster than some speed 
threshold, we divide the time by two. Conversely, we in¬ 
crease the time step by 20% if no expansion was needed 
during the last few time steps. The maximum time step 
is limited by the rate of change of the field and the prop¬ 
agation scheme, with the actual time step determined by 
the desired target accuracy. 

Note that the new algorithm can describe tunnelling 
as long as the accuracy is high enough to include the 
exponentially decreasing wavepacket inside the classically 
forbidden area. 


7. Implementation for multidimensional systems 

Consider a Hamiltonian of a system with spa¬ 
tial dimensions that is discretized on a multidimensional 
Fourier grid: 

Nd 2 

H 1 (g)... (g) (g)... (g) 1 -I- V{xi,...,xn^) , (7) 

ZlTli 
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where pi are N hy N matrices and Id is a by 
matrix, with N the number of points on the spatial grid 
of one dimension. For simplicity assume that the grid 
has the same number of points in each dimension. 

The kinetic part is separable; therefore each 1-d pi 
can be written analytically in the Fourier representa¬ 
tion, , then transformed to the unreduced B basis: 
pf = B~^piB. Here, as in Section Hi A 21 B denotes the 
matrix obtained via eq. [2] from the matrix G that con¬ 
tains the 1-d discretized Gaussians on the Fourier grid. 

Generally, the non-tensor product parts of the Hamil¬ 
tonian are the particle-particle interaction terms. We 
therefore decompose these particle-particle interactions 
as 

V (xi,X2 . . .Xat) 

mi TTln 

^ ■■■ (i) Vjl (a^i) ® {xn) 

tl=l jN=l 

( 8 ) 

where the single-particle potential vectors, are nor¬ 
malized, and c representing the magnitude of each term. 
The general problem of finding an optimal decomposition 
(with minimal possible error for any number of terms) is 
an open problem. For the purposes of this work we uti¬ 
lize the POTFIT algorithm which is optimal for 

two degrees of freedom and applicable generally. Each of 
the terms i^i) is then converted to the B basis. 

To obtain the matrix WB in eq. El one first defines 
the one dimensional set of discretized Gaussians, which 
gives a matrix G, then generates the corresponding one 
dimensional B matrix by i? = (G^)“^ and computes the 
elements of the two dimensional B matrix defined by the 
reduced basis btj = bi® bj. 

At this point, we have the whole Hamiltonian in the B 
basis, but stored in memory as a collection of matrices 
representing 1-d functions. The reduction process of the 
multidimensional B basis is then defined on the indices 
of the complete basis, icb S [1; One needs to store 

in memory a vector which maps the indices of the com¬ 
plete basis to the collection of 1-d indices of the reduced 
basis, in order to efficiently implement the many changes 
of reduced basis that occur during the computation of 
eigenmodes or the dynamics. 


8. Symmetry and Parallelizability 

Symmetry considerations can cut the number of 
Hamiltonian elements computed by two orders of mag¬ 
nitude [13 ■ Indeed, the B basis is a Gabor basis be¬ 
cause the G is Gabor by construction [43, which im¬ 
plies that it possesses a translation symmetry in phase 
space. From the latter, one can deduce that the in¬ 
tegral J bx^p^{xi,X 2 )*h{xi,X 2 )bx„p^ixi,X 2 ) depends on 
(pa — Pb) instead of (pa,Pb)- The exchange symmetry 


and the hermiticity of the Hamiltonian also reduce the 
number of elements to compute. 

Next, let us remark that PvB is straightforward to 
parallelize. Indeed, a significant part of the computa¬ 
tional effort is spent on converting elements of the re¬ 
duced Hamiltonian to the PvB basis — a task which is 
trivially parallelizable, as each element’s computation is 
independent. In practice, each time the reduced basis is 
expanded, only the new elements of the reduced Hamil¬ 
tonian need to be computed. For example computing 
the ground state on a grid x G [—100 a.u.; 100 a.u.] and 
|Ar| < 15 for an accuracy threshold W = 10“^ with six 
cores, we observe that over 90% of the task is linearly 
parallelizable. No significant resources have been spent 
to distribute the tasks between the cores or to merge back 
the results, which implies that this ratio will not decrease 
when the number of cores increases. This fact remains 
valid for a multi-machine architecture. The proportion 
of workload that is parallelizable is higher for dynam¬ 
ics than for the eigenmode problem, and increases for 
higher target accuracies, where the reduced Hamiltonian 
becomes larger. This can be compared, for example, with 
MCTDH that reaches only around 50% of parallelization 
for the dynamics on this type of grid because of the con¬ 
tributions from the SPF propagation that is not suited 
for the parallel MGTDH, as stated in [^ . 


9. Filtering of momenta correlations 

Ionization, produced by short NIR or XUV pulses, gen¬ 
erates travelling electron wavepackets. The momenta 
of these packets are actually lower than the maximum 
momentum component of the ground state wavepacket. 
Consequently, to analyze momenta correlations one has 
to remove the bound part of the electronic wavefunction, 
that would otherwise obfuscate the momenta correlations 
of ionized wavepackets. 

Within the phase space framework, this filtering op¬ 
eration becomes straightforward. Indeed, one just needs 
to define the phase space volume V corresponding to the 
bound states and remove from the reduced state the cor¬ 
responding 6-functions such that their orthogonal Gaus¬ 
sians are centered inside this phase space volume. The 
projector producing the filtering is: 

pf = 12 

feev 

where k G V means that the center of the fc-th Gaussian 
belongs to V. 

Here, to obtain the full momenta correlations plot, we 
define V by a simple spatial consideration, cutting every¬ 
thing with \xi\ < 30 a.u. However, a more subtle phase 
space definition of V could be use to analyze separate 
parts of the wavefunction. Phase space filtering includes 
coordinate space filtering as a special case. For example, 
one may filter out all single ionization components by 
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removing 6-functions such that the corresponding Gaus- 
sians satisfy the condition (|a;i| < Xr) or (|x 2 | < Xr), 
with Xr the radius of the atoms. This method does not 
require any additional computation, contrary to the one 
used in [l6| where many eigenstates need to be computed 
and then subtracted. 


B. Review of MCTDH 


and 


Then, following we define the single-hole function 




(fc) 




. /■* 


^Jk + l ^3N 


the mean-fields 


1. Sum-of-Products Ansatz 

Multi-Configuration Time Dependent Hartree 
(MCTDH, [^) is one of the leading approaches to 
solving the time-dependent Schrodinger equation with 
multiple degrees of freedom. The basic ansatz is as 
follows. Any multi-particle state in a finite Hilbert 
space may be decomposed as a sum-of-products of single 
particle functions (SPFs), 

if {XI,X2, ■ . .,XNi,t) 


mi m. 


= E- 

• E “A .-iNd 


{xNa,t) 


iiv=i 






(10) 

with 



II 


(x 2 ,t) I ^ = 0. Each term of this 

series is called a configuration. Under weak coupling 
conditions, the particles will be weakly correlated, and 
only a few configurations in the above expression will 
have non-negligible coefficients, allowing for a very effi¬ 
cient representation of the multi-particle wavefunction. 
For larger systems, it is beneficial for the decomposition 
to be made hierarchically, to match the hierarchy of 
system coinings, leading to a multi-layer MCTDH 
algorithm [^. For the purposes of this paper we made 
use of the original single layer MCTDH implementation 
package. 

To efficiently implement the MCTDH dynamics, one 
needs to decompose the Hamiltonian as well as the wave- 
function into a sum-of-products form. The decomposi¬ 
tion of the Hamiltonian is needed by both MCTDH and 
PvB, as it allows replacing multi-dimensional integrals 
with a series of one-dimensional integrals. Such integra¬ 
tions appear in both MCTDH and PvB when converting 
the Hamiltonian to the basis of interest. 


{H)f = (vE-f iFflTp)) 

and the density matrix 



Utilizing the Dirac-Frenkel variational principle, one ar¬ 
rives directly at 


L 







with the matrix of mean-fields. 

In MCTDH, eigenmodes are generally calculated by 
imaginary-time propagation, i.e. a relaxation technique. 
This algorithm has been extended to return multiple 
eigenmodes in a single run. 


C. Choice of pseudo-spectral basis 

Both PvB and MCTDH require an underlying dis¬ 
crete basis of localized functions to represent the wave- 
function at the grid points. More precisely, MCTDH 
uses exclusively one-dimensional grids that support the 
SPFs whereas the PvB lattice covers the multidimen¬ 
sional phase-space but is never entirely used. While many 
possible pseudo-spectral bases are pos sible for MCTDH 
and PvB, as shown in [s^ and [44| respectively, we 
have chosen the Fourier pseudospectral functions for PvB 
and the harmonic oscillator pseudospectral functions for 
MCTDH for the sake of simplicity. 


D. Scaling of PvB vs. MCTDH 


2. Dynamics 

MCTDH dynamics requires equations of motion for 
both the tensor as well as the SPFs, Let us 

start by introducing the following simplified notations: 


Let us compare the numerical effort required for 
PvB vs. MCTDH. With MCTDH, the effort has two 
terms, one from the basis function evolution and one 
from the calculation of the coefficients: Nt{mndN'^ + 
where Nt is the number of time steps, m 
is the number of terms in the POTFIT series, n is the 
number of single particle functions in each dimension, d 
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is the dimension and N is the number of points in the 
one dimensional grid. For low dimensional problems on 
large grids the first term NtmndN'^ dominates. With 
PvB the numerical effort decomposes into two parts. The 
first is the cost to precompute the Hamiltonian terms in 
the PvB basis, mNrN'^ + NrN^, where Nr is the num¬ 
ber of functions in the reduced basis. The first term 
mNrN"^ corresponds to the conversion of the POTFIT 
series and the second term NrN^ comes from the conver¬ 
sion of the kinetic Hamiltonian. The second part of the 
effort is the dynamics itself, scaling as NtN^. Note that 
since a large number of timesteps Nt is required to ade¬ 
quately represent the control field, the term proportional 
to Nt dominates for the calculations in this paper. For in¬ 
stance in the example used in Section HIC, Nt = 15000 
and 5000 < Nr < 10000 while N = 4000, n = 12 and 
m K, N. Considering just the term proportional to W, 
this translates to a two order of magnitude speedup of 
PvB as compared to MCTDH, including the cost of the 
conversion to the PvB basis. In the case where the con¬ 
version results are already available from a previous run, 
the speedup increases to four orders of magnitude. This 
difference in speed between the two methods increases 
with the size of the grid. 



FIG. 2. POTFIT: Sum-of-products representation for the 
Coulomb potential Vapprox on a grid x € [—100,100] with 
300 points compared to the exact potential V. Note that the 
removal of even a single term from the expansion represent¬ 
ing the Coulomb decomposition leads to a signihcant drop in 
accuracy. 

The Hamiltonian used is therefore: 


Note that if the POTFIT series were shorter, MCTDH 
dynamics would scale very similarly to PvB dynamics. 
Thus, this large advantage of PvB is valid only for po¬ 
tentials that are difficult to decompose, which is the case 
for the two-electron Coulomb potential as explained be¬ 
low. 


III. APPLICATION TO 1-D HELIUM 


H = 


47reo 


-b 


QeQ 


y^xl+a, 
^2 


+ 


+' 


\J[xi - X2f + aj 


-b 


2m. 


{pI+pI) (11) 


with Q = —2qe, qe being the electron charge. Note that 

= , is notoriously diffi- 


the interaction term, , ^ —=—=, 

cult to represent as a sum-of-products, requiring a very 
large number of elements in the series to achieve an ac¬ 
curate representation (see fig. [2]). 


A. The Model System 


B. Ground State 


To test the suitability of PvB and MCTDH for simula¬ 
tions that combine both bound and unbound dynamics, 
we apply both methods to the double ionization of ID he¬ 
lium. We shall begin with the calculation of eigenmodes, 
and then proceed to the ionization dynamics. 


Our benchmark system, the ID helium model, consists 
of two electrons, each with a single degree of freedom, in¬ 
teracting with each other and a central (nuclear) poten¬ 
tial. We use a regularized form of the Coulomb potential, 
= , where the regularizer, oq, is set to 0.739707902, 


such that the ground-state energy of the model matches 
the experimen tally measured binding energy of helium, 
2.903385 amu [^, [^. For the purpose of high accuracy 
benchmarking, some of the following results take as ref¬ 
erence a ground state energy E = —2.90338599 where 
only the first six decimals are experimentally relevant. 


1. Computational Speed and Memory 

We begin by comparing the relative efficiency and be¬ 
havior of MCTDH and PvB when solving the time in¬ 
dependent Schrddinger equation. In the following dis¬ 
cussion the range x G [—100 a.u. ; 100 a.u.] is used for 
calculating the ground-state. 

We used the so-called improved relaxation method im¬ 
plemented in MCTDH to solve the eigenvalue problem. 
The MCTDH accuracy improves rapidly with the num¬ 
ber of configurations, as shown in fig. |3] where 80 con¬ 
figurations are enough to reach 10“® accuracy for the 
ground energy. Note that we consider here the total 
number of configurations and not the number per di¬ 
mension. PvB converges more slowly as a function of 
the basis size: for the same accuracy 3500 cells are re¬ 
quired, as can be seen in fig. S) Thus MCTDH is gen¬ 
erally faster than the phase-space algorithm by an order 

















FIG. 3. Convergence of the MCTDH ground state energy 
as a function of number of configurations. 



FIG. 4. Convergence of the PvB ground state energy as 
a function of reduced basis size. 


ergy oi E = —2.903379690969, where the improved re¬ 
laxation converged up to the last digit of this value. On 
the other hand, the PvB algorithm halts with an accu¬ 
racy estimate of 10“^, because it detects that the ground 
state has reached the edge of the phase-space area. This 
accuracy is in fact the real accuracy, as can be verified by 
running the MCTDH improved relaxation with N = 255 
points, which produces E = —2.90338601. 

Insight into the problem can be obtained by inspect¬ 
ing the phase-space representation of the state. Figure [5] 
shows the one dimensional projection (partial trace fol¬ 
lowed by summation of absolute value amplitudes) of 
the MCTDH final state represented in the PvB basis. 
Each unit square represents the position of a Gaussian in 
phase-space and the color map represents the amplitude 
of the overlap between this Gaussian and the state of 
interest. The upper and lower plots correspond respec¬ 
tively to = 100 and N = 255. The red lines depict 
the limits of the phase-space in the case N = 100. In the 
case N = 255 the state wavenumber goes well beyond the 
N = 100 limit, which is another way to say that N = 100 
has insufficient resolution to describe the high frequency 
components of the ground state, leading to inaccurate 
results. For a large grid, e.g for ionization problems, the 
long running times make it important to avoid repeating 
the same computation with different grid sizes to ensure 
convergence. 


of magnitude. However it requires a tuning of the initial 
guess, whereas the phase-space algorithm only needs to 
choose the grid size and the desired accuracy. Moreover, 
it is interesting to compare how many complex numbers 
are needed to store the state in memory. In the case of 
MCTDH, it corresponds to the number of configuration 
times the size of the one dimensional grid, which gives 
here 80 x 300 = 24000, whereas for PvB it is given di¬ 
rectly by the reduced basis size: 3500 in this case. There¬ 
fore, for low dimensional problems, the representation of 
the state is more efficient in the PvB method. 


2. Accuracy Control 

One of the distinct advantages of the phase-space al¬ 
gorithm is its accuracy control. This is in contrast to 
MCTDH, which is unable to detect when the number 
of grid points is insufficient. Several runs with different 
grid resolutions are necessary to ensure convergence. If 
the number of points is too small, MCTDH’s improved 
relaxation will converge with machine accuracy towards a 
wrong eigenvalue. In contrast, the phase-space algorithm 
knows intrinsically which accuracy it achieves. 

To illustrate this point, consider a grid with N = 100 
points in each dimension, with x G [—20 a.u. ; 20 a.u.]. 
This corresponds to a maximal k of fc^ax = Nn/L = 
7.85. On this grid, MCTDH returns a ground-state en- 



FIG. 5. Comparison of the one dimensional phase-space pro¬ 
jections of the ground state, obtained with N = 100 grid 
points (upper) and N = 255 (lower). 


The PvB algorithm gives not only an a posteriori evalu¬ 
ation of the accuracy; it allows an a priori estimate. This 
estimate depends only on the Wavefunction Amplitude 
Cutoff, W, which determines which cells are included in 
the reduced basis, based on the value of \{gk\ ■*/')!■ Con- 
















9 


sider the reduced state and its complement 


Nr _ N 

'ip'j =^Ck\bk) and •0^= ^ Ck\gk) ■ (12) 

k—1 k—NR-\-l 


The state then decomposes to \tjj) 


ill 


Ip). Assuming 


|^/>) is an eigenstate of energy E, we can quantify the 
contribution of the two sets to the energy: 


E={P;\H\'iP) 

= E^°'> + 2e(^^ (P) H P) 

Nn, 

= ^ c*Ck {2E {g,\ gk) - {g,\H \gk)) (13) 


where E^^'^ = (pj 


H p)j and is the number of ele¬ 
ments in the complement set that are close enough to 
the border of the reduced set to have non negligible 
coefficients. _We have also made use of the fact that 
Ip tp\ = (ip pj), as the reduced set and its complement 


are orthogonal. Note that there is no first order term in 
Cfe, indicating that a perturbation to the eigenstate does 
not contribute to the energy to first order. This is as 
expected from the Rayleigh-Ritz variational principle. 

By construction, the coefficients in the second term 
satisfy Ck < W. Thus, the energy error is bounded by 
a term of the form \E — where the 

constant C bounds the term {2E {gj \ gk) — {gj \ H l^fc)). 

From these considerations, we can deduce a heuristic 
bound on the error for the ID helium model. We first 
note that Nyj itself depends on W because a larger re¬ 
duced basis is needed to run a lower W computation. For 
the ID helium model we observe in the numerical simula¬ 
tion that Ntj, scales approximately as Thus, the 

overall bound error bound on | FI — scales approxi¬ 
mately as as shown in fig. [SI A systematic study 

would be required to validate this bound for others sys¬ 
tems, however W itself can still be used as a conservative 
estimate of the error in other systems. 


C. Dynamics 

We now proceed to the dynamics calculations, per¬ 
forming all calculations in the velocity gauge. We first 
validate the PvB code by comparing to MCTDH on a 
small grid. Figured shows the electronic wavefunction at 
different times under the influence of just an NIR pulse. 
In the velocity gauge, the controlled Hamiltonian takes 
the form: 


H = Hq —+ ux\jv{t)){pi + P 2 ) , (14) 

TOe 

where Hq is the drift Hamiltonian (eq. fTTl) and u is 
the amplitude of the electric potential. For this first 



FIG. 6. Error in the ground state energy as a function of the 
wavefunction amplitude cutoff W. 


test, the XUV pulse is set to zero. The NIR pulse 
is taken to have a sine envelop in order to have ex¬ 
actly zero derivatives at the beginning and at the end: 
MNiR = ANiRsin(27rt/rNiR - 7r)sin(7rt/(4rNiR))^ with 
Tnir = 110.32 a.u. (= 2.6685 fs). This corresponds to 
a wavelength of 800nm, and a total duration of 10.67fs. 
The peak amplitude is Anir = 0.6627 (corresponding to 
an intensity of 5 x lO^^W/cm^). 

The computation is carried out for an x range 
[—100 ... 100] and N = 1000 points in each dimension. 
Note that what is shown is a one dimensional projection 
on the Fourier grid of the two dimensional state. Al¬ 
though the PvB state is pixelated, the projection onto 
the Fourier grid is smooth, since the PvB state contains 
by construction exactly the same information. 

In order to obtain a rich ionization dynamics, we add 
to the NIR pulse two XUV pulses. While it has long been 
known that helium can be ionized by a strong NIR pulse, 
the combination of an NIR pulse to excite the atom to 
high bound states, in conjunction with one or more XUV 
pulse triggering the ionization, allows for much better 
control of the resulting dynamics. We based these pulses 
on the one used in [l6| which has the form: uxuv = 
Axuvsin(27rt/rxuv)exp(-(t - 5rxuv/4)^/(2(T^)) with 
2xuv = 2.07 a.u. which corresponds to a wavelength 
of 15 nm. We take cr = 6.207 a.u. (150 attoseconds) and 
peak amplitude Axuv = 0.00176 a.u. (corresponding to 
an intensity of = 1 x lO^^W/cm^). However, to obtain 
a large amount of double ionization and make the quali¬ 
tative analysis clearer we increased the amplitude of the 
two XUV pulses by a factor of 50. The complete control 
pulse is shown in fig. [HI The control strategy is to gener¬ 
ate a sequential double ionization with the two successive 
XUV pulses. For this simulation the PvB calculation is 
based on an underlying grid with a range x G [—400; 400] 
with N = 4000 points in each dimension. On such a 
grid, the dimension of the unreduced Hilbert space is 
16 X 10® while the maximum dimension of the reduced 
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X Time (a.u.) 


FIG. 7. Comparison of the one dimensional projection of 
the dynamic state computed with PvB and MCTDH under 
the action of the NIR pulse at two different times (top and 
bottom rows). Note how the center peak of the wavefunction 
is tilted left and right, depending on the phase of the driving 
field, (the red dot in the right panel). Also note the excellent 
match between the two methods. 


Hilbert space used during the dynamics is 28207. This 
translates to a reduction by 5 orders of magnitude in the 
number of elements in the Hamiltonian, as compared to 
the size of the unreduced Hamiltonian. 


With these settings, MCTDH is typically two orders 
of magnitude slower than PvB. The difference in speed 
is a result of the large number of terms required in the 
POTFIT series when dealing with the long range multi¬ 
electron Coulomb potential. For this range, the number 
of terms in the POTFIT series becomes on the order of 
N (see fig. ED, which drastically reduces the efficiency of 
MCTDH, as explained in subsection III PI 


The two electron dynamics are shown in figs. [3 [9] and 
[TUI Figures [3 and [TU] have the following structure. Frame 
(a) is the wavefunction projected onto the {xi,X 2 ) plane. 
Single ionization corresponds to the wavepacket moving 
along the horizontal and vertical axes of the frame (a) 
while double ionization corresponds to the wavepacket 
moving along the diagonals of the same frame. The 
frames (b) and (c) show the phase space projection of 
the two electrons. These frames are identical due to sym¬ 
metry, up to a switching of axes to facilitate the reading 
of the correlations. Frame (d) presents the wavefunc¬ 
tion projected onto the {pi,P 2 ) plane, after being fil¬ 
tered to remove the contribution from the bound states. 
The plots are produced by a 2-d projection of the 4-d 
phase space amplitude, e.g. the amplitude of one cell 

in the spatial correla¬ 


1 ^ 1 = (^gix[''\p["'\x''2\p^^'^)\tp 


tions plane is: 


C 




\ 


E 




(15) 


where {x^^\p^ 4 ^\x 2 \p^^) denotes the centers of the 
Gaussians. This choice may not be the best as it induces 
a strong loss of information but is is straightforward to 
implement numerically and efficient enough for the pur¬ 
pose of analyzing the position of the wavepackets in phase 
space. 

Figure [3 shows the projections at t = 232 a.u., be¬ 
fore the second XUV pulse. Inspection of the wavepacket 
dynamics reveals that the amplitudes marked by B and 
C are generated by the NIR pulse while the amplitude 
labelled by D has been generated by the first XUV 
pulse. This labelling is confirmed by checking the mo¬ 
mentum expected from the frequency of the pulses with 
the de Broglie relation, which gives pxuv ~ 3.04 and 
PNiR ~ 0.05. Of course, one should not expect to ob¬ 
serve exactly these values, since the dynamical process 
is complicated, but they allow one to discriminate be¬ 
tween travelling wavepackets with a low |p| and higher 
IpI, produced respectively by the NIR and XUV pulses. 
Note that the projections (dashed blue lines) in frames 
(b) and (c) allows one to distinguish the different contri¬ 
butions to the double ionization wavepacket. The pro¬ 
jection lines show that the double ionization amplitude 
may include contributions from both pulses. Frame (d) 
shows the expected momenta of the different wavepack¬ 
ets. Note that after the bound amplitude is filtered out 
all momenta are lower than the maximum momentum of 
the bound states, that are not filtered out in frame (b) 
and (c). 

Figure [3 shows three snapshots to emphasize the pres¬ 
ence of concerted and sequential double ionization. The 
middle column presents the same projection as frame (c) 
of figures 15] and [TUI The left column plots the {xi — X 2 ) 
correlations, while the right column plots the (pi — P 2 ) 
correlations filtered such that only the double ioniza¬ 
tion is visible. On the first line, only single ionization is 
present, thus the right plot on this line is empty. The sec¬ 
ond line shows the first appearance of double ionization 
{t = 219.7), which corresponds to only one cell in the 
first column, indicated by the green circle at xi ~ X 2 - 
One can see in the right column that both momenta 
are nonzero, and approximatively equal. The combina¬ 
tion xi ~ X 2 and Pi « p 2 is the signature of concerted 
double ionization. Thus, there is a small amplitude of 
concerted ionization at this point in time. The double 
ionization becomes greater at later time and largely se¬ 
quential, as shown in the third line (t = 274.7). The 
green rectangle in the left frame of the third line is lo¬ 
cated at 11x211 ||a:i||. The right frame of the third line 

shows the filtered momentum correlations corresponding 
to that rectangle. We observe the wavepacket, originat¬ 
ing at xi Ri 0, X 2 < 0, has momentum p 2 ~ 0 with a 
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FIG. 8. Two electron wavefunction before the second XUV pulse at t = 232 (a.u.). A double ionization marked by (A) is 
generated by the influences of NIR pulse, point (B) and point (C), and the contribution of the wavepacket ionized by the first 
XUV, marked by (D) 


significant negative pi. This is the signature of sequen¬ 
tial double ionization. The green arrows on the first col¬ 
umn show the direction of the motion obtained from the 
average of the filtered wavepackets shown on the third 
column. 

Figure [TOl shows the projections at i = 275 a.u., fol¬ 
lowing the second XUV pulse. The second XUV pulse 
generates a significant amount of sequential double ion¬ 
ization (fig. [TOk . points A and C) as well as a second 
wave of single ionization (fig. [TPk . point B). The double 
ionization in this case is sequential, because the corre¬ 
sponding wavepackets arrived at their present location 
by a two step process, first following the axis before leav¬ 
ing it to get closer to the diagonal and anti-diagonal of 
fig. [TOb . The NIR pulse has the effect of shifting the new 
single ionization wavepacket toward positive momentum 
compared with the single ionization from the first XUV. 
Indeed, it can be seen in fig. [THh that the two lobes with 


higher \x 2 \ generated by the first XUV pulse have lower 
momentum than the one generated by the second XUV 
pulse. Note again that the dashed blue projection lines 
allow one to determine by eye the contributions to the 
double ionization amplitude. For example the amplitude 
marked C corresponds to a superposition of amplitude 
from G, generated by the first XUV pulse, from E, gen¬ 
erated by the NIR pulse and from D, generated by the 
second XUV pulse. 

Figure [TT] shows the photoelectron momentum distri¬ 
bution at the end of the pulse. We observe the expected 
contributions from the different parts of the pulse. The 
correlations are bounded by the momenta of highest ab¬ 
solute value, corresponding to the XUV pulse. For a 
given electron, each of these high momentum wavepack¬ 
ets correlates with all momenta of the other electron, 
which creates the square shape. One can also notice the 
ponderomotive shift between the two XUV wavepackets. 
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FIG. 9. Snapshots of the two electron wavefunction between the two XUV pulses, at times t = 215.6 a.u., t = 219.7 a.u. 
and t = 274.7 a.u., from top to bottom. The first column shows the xi — X 2 correlations. The second column shows the one 
dimensional phase space projection of the first electron. The third column shows the filtered pi — p 2 correlations, see text for 
details. The times corresponding to the three snapshots are indicated by red circles on the control pulse in the bottom frame. 
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FIG. 10. Two electron wavefunction after the second XUV pulse at t = 275 (a.u.). A sequential double ionization marked by 
(A) and (C) is generated by the influence of the second XUV pulse marked by (D) and the NIR marked by (E), which correlates 
with the first ionization marked by (F) and (G). The second XUV also generates another single ionization wavepacket, marked 
by (B). 


for example for P 2 = 0 and pi > 0. This shift is approx- 
imatively 2 atomic units, which corresponds to the NIR 
total amplitude. Note that the drop in intensity of the 
boundary layer comes from the phase space truncation, 
these cells having some contributions from the core cells 
that were filtered out. 

Finally, note that the accuracy of this simulation is 
10“® for li/jp, {W = 10“'’’), which is similar to many 
other studies, see for example m- 


IV. CONCLUSION 

We have applied the PvB method to the calculation of 
both eigenstates and ionization dynamics. The ground 
state calculation (and eigenstate calculations in general) 
requires a relatively small grid, and in such cases PvB is 


slower than MCTDH. However, we have shown that PvB 
offers better control of the accuracy. On the dynamics 
side, we showed that PvB and MCTDH reproduce the 
same short-time dynamics, but for the large and dense 
grids required for the long-time dynamics PvB is faster 
than the Heidelberg MCTDH code by orders of mag¬ 
nitude. Indeed, it can deal more easily with the high 
precision propagation required to simulate the low am¬ 
plitude wavepacket that leaves the nucleus. Moreover, 
the phase space representation allows one to identify by 
eye the mechanisms leading to the formation of the dif¬ 
ferent correlated wavepackets. Finally, PvB computa¬ 
tions are straightforward to parallelize to multi-cores and 
multi-machine environments. PvB also requires very lit¬ 
tle initial tuning: only the desired accuracy and the phase 
space size have to be defined. 

The main current limitation of PvB lies in the size 
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FIG. 11. Photoelectron momentum distribution at the end of the pulse. 


of the reduced basis. This limitation can potentially be 
overcome by decomposing the basis in terms of sums of 
product of one dimensional terms. Therefore, it is inter¬ 
esting to consider the possibility of merging PvB with 
MCTDH. We note that the two methods, MCTDH and 
PvB, use fundamentally different approaches: MCTDH 
aims at separating the degrees of freedom while PvB fo¬ 
cuses on reducing the effective grid size. Since the two 
methods do not have the same domains of efficiency, it is 
very possible that merging the two methods could pro¬ 
vide a method with the advantages of both approaches, 
able to simulate both high dimensional and long range 
dynamics. For instance, one may consider replacing 
the one dimensional discrete variable representation of 
MCTDH by the PvB method, or alternatively, using PvB 
to describe each electron in 3-d and introducing the cor¬ 
relation via MCTDHF, i.e. MCTDH with exchange sym¬ 


metry built in via the Slater determinants. 

The PvB code developed for this project is available 
upon request and we would be happy to cooperate in 
further developing the methodology and testing new ap¬ 
plications. 
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